% This code calculates the subjective risk premia using a Gaussian 
% quadrature approach.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f = calculate_subRet_proj(d_c_diff, tmu, beta_PD_proj, I, v, fspace_dc)

% INPUTS: 
% -- d_c_diff: d_t - c_t, one of the state variables
% -- tmu: \tilde{mu}_t, one of the state variables 
% -- beta_PD_proj: Coefficients of polynomials from P/D approximation
% -- I: Parameter structure for other parameters
% -- v: Number of points for intergral
% -- fspace_dc: Function space that contains all the projection
%               information for d_t - c_t
% OUTPUT:
% -- f: Subjective returns as in Equation (D.40)


%% Some coefficients for use later
coeff_1 = sqrt(1 + I.nu) * I.sigma;
coeff_2 = coeff_1 * I.nu;
coeff_3 = I.lambda * sqrt(1 + I.nu) * I.sigma;

%% Gaussian quadrature sampling. Taking +/_ 8 s.d.
[epsilon_tmu , weight_tmu] = qnwlege(v,-8,8);             %For epsilon from -8 to 8, do Gaussian-Quad sampling
[epsilon_dc , weight_dc] = qnwlege(v,-8,8);               %For epsilon from -8 to 8, do Gaussian-Quad sampling

%% Calculate subjective returns
%--------------------------------------------------------------------------
% For the following calculations see Appendix D.4 for formulas
%--------------------------------------------------------------------------
% Calculate the H matrix, which is the log P/D ratio
grid_H{1,1} = tmu;
grid_H{1,2} = d_c_diff;
exp_H_matrix = exp(funbas(fspace_dc, grid_H)*beta_PD_proj) * ones(1 , v^2);

% Next-period log P/D ratio
grid_H_new{1,1} = tmu + coeff_2 * epsilon_tmu;
grid_H_new{1,2} = (1 - I.alpha) * d_c_diff + I.alpha * I.mu_dc + (I.lambda - 1) * (tmu + coeff_1 * epsilon_tmu) + I.sigma_d * epsilon_dc;
H_matrix_new = (funbas(fspace_dc, grid_H_new)*beta_PD_proj)';
exp_H_matrix_new = exp(H_matrix_new);

% Shocks to the 2-D normal distribution
Shock_density_matrix = (0.5/pi) * exp(-0.5 * (kron(ones(v , 1) , epsilon_tmu).^2 + kron(epsilon_dc , ones(v,1)).^2)');

% Additional shocks from the dividend growth term
Add_Shock_Matrix = exp((kron(ones(v,1) , coeff_3 * epsilon_tmu) + kron(I.sigma_d * epsilon_dc , ones(v,1))))';

% Multiply everything together
Expectation_integral = (Shock_density_matrix .* Add_Shock_Matrix .*((exp_H_matrix_new + 1)./exp_H_matrix)) * kron(weight_dc , weight_tmu);

% Finally multiply the constant terms from dividend growth
f = exp( I.lambda * tmu - I.alpha * (d_c_diff - I.mu_dc) ) * Expectation_integral;
end